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We consider a Bose-Einstein condensate (BEC), which is characterized by long-range and 
anisotropic dipole-dipole interactions and vanishing s-wave scattering length, in a double-well po- 
tential. The properties of this system are investigated as functions of the height of the barrier 
that splits the harmonic trap into two halves, the number of particles (or dipole-dipole strength) 
and the aspect ratio A, which is denned as the ratio between the axial and longitudinal trapping 
frequencies ui z and uo p . The phase diagram is determined by analyzing the stationary mean-field 
solutions. Three distinct regions are found: a region where the energetically lowest lying stationary 
solution is symmetric, a region where the energetically lowest lying stationary solution is located 
asymmetrically in one of the wells, and a region where the system is mechanically unstable. For 
sufficiently large aspect ratio A and sufficiently high barrier height, the system consists of two con- 
nected quasi-two-dimensional sheets with density profiles whose maxima are located either at p = 
or away from p — 0. The stability of the stationary solutions is investigated by analyzing the Bo- 
;y>} 1 , goliubov de Gennes excitation spectrum and the dynamical response to small perturbations. These 

studies reveal unique oscillation frequencies and distinct collapse mechanisms. The results derived 
bJ)' within the mean-field framework are complemented by an analysis based on a two-mode model. 
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Dipolar BECs have recently attracted a lot of attention both theoretically and experimentally [j], Q • The experi- 
mental realization of a Cr BEC just a few years ago constitutes an important milestone 3]. Compared to alkali atoms, 
Cr has a comparatively large magnetic dipolc moment of 6/is , which leads to an enhancement of the dipole-dipole 
interactions by a factor of 36 compared to alkali atoms. The anisotropy of the dipole-dipole interactions has been 
observed experimentally by analyzing time of flight expansion images of Cr BECs released from a cylindrically sym- 
metric external confining potential [j]. If combined with theoretical calculations, the time of flight images reveal the 
initial density distribution of the dipolar gas and depend, e.g., on whether the magnetic dipole moments are aligned 
along the axial or longitudinal confining directions, respectively. Furthermore, by taking advantage of the tunability 
of the s-wave scattering length near a magnetic Fano-Feshbach resonance, the relative importance of the dipole-dipole 
interactions can be changed [1,(1,01, paving the way for a variety of interesting experimental studies. Dipolar BECs 
are characterized by intriguing collapse mechanisms 1, H, H, H, fiol . ITTI . [l2l [l3l Il4l . uBi . [H, [ItJ , and unique excitation 
spectra [H, [H, and vortex structures [20 L 21. 22 ri23l |24[] . In addition, dipolar gases loaded into optical lattices 
may allow for the realization of novel phases [25l . 

Although the dipole-dipole interactions in alkali gases are too weak to result in observable effects in most experi- 
ments, it is believed that they play a decisive role in the formation of spin textures in 87 Rb condensates [H, [29|, in 
the dynamics of Bloch oscillations of 39 K BECs loaded into an optical lattice [3C| and in 7 Li BECs 1311 . Furthermore, 
BECs and degenerate Fermi gases that consist of polar molecules may be realized in the near future|32l. [33. 34]. This 
prospect adds a new intriguing twist since the interactions between two polar molecules can be tuned by an external 
electric field [11] . This opens the possibility to enter the strongly-correlated regime and thus to realize a variety of 
condensed matter analogs [H, [H, HtJ • 

This paper considers an aligned dipolar BEC in a double well geometry. Double well potentials play an important 
role in chemical and condensed matter physics, among other areas. In the context of cold atom physics, s-wave 
dominated alkali systems in a double- well have, e.g., been used to study Josephson-type oscillations (3a. [33. l4fl l4ll. 
H, EH, 0, |H, 111] . The density oscillations of the Bose gas can be interpreted as corresponding to the charge current 
that characterizes "standard" condensed matter Josephson junctions. Related to this, the macroscopic quantum 
self-trapping of atoms in one of the wells has been demonstrated experimentally and has been interpreted within a 
two-mode model that can be derived from the Gross-Pitaevskii (GP) equation [44j, |45| . The double well system has 
also been used to experimentally study spin-squeezing [471 ]. In this context, the left and the right wells of the system 
serve as the two arms of an interferometer [48j . The number difference and relative phase of the double well system 
are conjugate variables, whose combined measurements has revealed that the system is entangled [471 ]. 

Here, we investigate the behaviors of aligned dipoles under cylindrically symmetric harmonic confinement with 
a repulsive Gaussian potential centered at z = 0. We limit ourselves to situations where the dipoles are aligned 
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along one of the symmetry axis of the external harmonic confining potential. This restriction reduces the parameter 
space and also significantly reduces the numerical efforts. Arguably, it may be the conceptually simplest case. We 
are particularly interested in determining the phase or stability diagram for both cigar-shaped and pancake-shaped 
harmonic confinement. The boundaries of these phase diagrams are governed by the non-trivial interplay of the dipole- 
dipole interactions, the energy due to the external harmonic confinement, the energy due to the Gaussian potential 
and the kinetic energy. The interplay of these energy contributions leads to density profiles unique to anisotropic 
interactions. In agreement with Ref. (49j . we observe Josephson oscillations as well as macroscopic quantum self- 
trapping of the system for appropriately chosen parameters. We characterize the transition between these two regimes 
by analyzing the excitation spectrum and the real time response of the system to a small perturbation. Our study 
of two neighboring pancake-shaped dipolar gases can be viewed as a first step towards understanding a multi-layer 
system of dipolar pancakes. For a single pancake, an angular roton instability has been predicted to occur [l5j]. For a 
layer of two-dimensional dipolar BECs, a new length scale is given by the interlayer distance and the roton instability 
has been predicted to be enhanced compared to the single layer case [Hfj. Other multi-layer studies can be found in 
Refs. [HI, El. 

The remainder of this paper is organized as follows. Section HT1 introduces the mean-field GP equation and discusses 
how we determine the stationary and time-dependent solutions. The Bogoliubov de Gennes equations that are 
employed to determine the excitation spectrum of the dipolar gas are introduced. Section Mil reviews the two- mode 
model which provides an intuitive understanding of the time-independent and time-dependent GP solutions in the 
small A regime. Section IIVI presents our stationary solutions. We discuss the phase diagram as functions of the 
number of particles (or equivalcntly, the mean- field strength), the aspect ratio and, in selected cases, the barrier 
height. Section [V] presents our time-dependent studies. We investigate certain dynamically stable regimes and 
deduce distinct collapse mechanisms from the response of the system to a small perturbation. In addition, selected 
Bogoliubov de Gennes eigenmodes are discussed. Lastly, Sec. IVII summarizes our main findings and discusses possible 
future studies. 



II. MEAN-FIELD DESCRIPTION OF DIPOLAR BECS 



Section III Al introduces the time-dependent mean- field GP equation for a dipolar BEC and discusses the numerical 
techniques employed to determine stationary and time-dependent solutions. Section III Bl reviews the Bogoliubov de 
Gennes equations for the dipolar BEC. 



A. Gross-Pitaevskii equation 

The time-dependent GP equation for a dipolar BEC consisting of N identical point dipoles is given by [H, [T3, HH 

i»^=2W,t), (1) 

where the mean-field Hamiltonian H reads 

H = -^V 2 + F cxt (r) + 
2m 

(N-l)f V dd (r- ^(r 1 , i)| W. (2) 

Here, m denotes the mass of the dipoles. We interpret ip(f,t) as a single-particle wave function and correspondingly 
use the normalization J \ip(r, t)\ 2 d 3 f = 1. The external cylindrically-symmetric confining potential V cyi t consists of a 
harmonic trapping potential \\\ with angular frequencies uj p and uj z and a Gaussian barrier V g with height A (A > 0) 
and width 6, 

F cxt (f) = F ho (p,2) + V g (z), (3) 

where 

V ho (p,z) = -m(uy + u 2 z z 2 ) (4) 

and 



F g (z)=Aexp(-— ). (5) 
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We define the aspect ratio A of the harmonic confining potential as 

A = £. (6) 

UJ p 

Throughout, we employ cylindrical coordinates and write r = (p, (f, z). 

The third term on the right hand side of Eq. @ represents the mean-field potential, which depends on both the 
density of the system and the dipole-dipole potential V^d- Throughout, we assume that the dipoles are aligned along 
the z-axis, 

_, N , 2 l-3cos 2 $ 
V Ad {r-r*) = d 2 7 
\r — r\ 6 

where d denotes the strength of the dipole moment of the dipolar atom or molecule under study and $ the angle 
between the relative distance vector r — r 1 and the z-axis. Throughout, we assume that the s-wave scattering length 
a s vanishes, implying the absense of the usual s-wave contact interaction term in Eq. (|2|). For dipolar Cr BECs, e.g., 
this can be achieved by varying an external magnetic field in the vicinity of a Fano-Feshbach resonance 0, 0] ■ 
Rewriting the integro-differential equation [Eq. fT} with Eqs. flU)-©] in harmonic oscillator units a z and E z , where 



moj z 
and 



(8) 



E z = Tiw z , (9) 

shows that the GP equation depends on four dimensionless parameters: i) d 2 (N — 1)/ (E z a z ), which characterizes 
the strength of the mean-field potential; ii) the aspect ratio A; iii) the scaled barrier height A/E z ; and iv) the scaled 
barrier width b/a z . To reduce the parameter space, we consider a fixed barrier width b, i.e., b — a z /5. While most of 
our calculations are performed for A = 12E Z , we consider smaller barrier heights in selected cases. The aspect ratio 
is varied from A = 0.1 (cigar-shaped external harmonic confinement) to A = 12 (pancake-shaped external harmonic 
confinement). Lastly, the dimensionless mean-field strength D, 

d 2 (N-l) , , 

E z a% 

is, for a given A and A, varied from to the value D cl - at which collapse occurs. 

In practice, the mean-field strength D can be adjusted by loading the double-well potential with condensates of 
varying particle number TV. More conveniently, one might envision tuning the electric dipole moment of a molecular 
BEC through the application of an external electric field [35[ or, in the case of magnetic Cr BECs, by changing the 
ratio between the dipole-dipole and the s-wave interactions through the application of an external magnetic field in 
the vicinity of a Fano-Feshbach resonance 0,S0|- Although our study considers a s — and varying D, the latter 
scenario should allow for the observation of a number of features predicted in this study. Experimentally, the Gaussian 
barrier potential of varying height and width can be realized by a repulsive dipole beam with adjustable intensity and 
waist. 

The solutions to the integro-differential mean-field equations have to be determined self-consistently since the 
density |^| 2 , which is part of the solution sought, also enters into the mean-field potential. The stationary solutions 
can be written as ip(r) — ^{p, z)h(ip) with h(ip) = exp(jfcy) / 'y/2n . In the following, we seek stationary solutions with 
azimuthal quantum number k — 0. Our calculation of the excitation spectrum do, however, include k > modes 
(see Sec. Ill Bp . The evaluation of the integral contained in the mean- field potential can be performed most readily 
by transforming to momentum space via a combined Fourier-Hankel transform (54[. To determine the stationary 
solutions, we implemented two different approaches: i) We minimize the total energy of the system following the 
conjugate gradient approach (55j. In this approach, the solution is expanded in terms of harmonic oscillator basis 
functions in p and z, and the expansion coefficients are optimized so as to minimize the total energy per particle, ii) 
We propagate an initial state in imaginary time till the stationary solution has been projected out. 

The basis functions and the initial state are both represented on a grid in the p and z directions. The grid along p 
is chosen according to the zeroes of the Bessel functions (see Ref. [Hj]), which are distributed roughly linearly. Along 
the z-direction, we use a linear grid. For most calculations, a grid of N p x N z = 64 x 128 is sufficient. We employ a 
rectangular simulation box of lengths [0,/9 max ] and [— z max , z max ]- For pancake-shaped systems (i.e., A > 1), a "cutoff" 
is used for the dipolar potential (i.e., the interaction is truncated for \z\ > z max ), which reduces the interaction of the 
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true BEC with an "artificial image BEC" and thus allows for the usage of a smaller z max [H[. F° r A < f , no cutoff is 
employed. Typical values for p max and z max are around 15a p and 12a z , respectively, where a p — J H/ {muj p ) . 

We have checked that the conjugate gradient and imaginary time evolution approaches result, within our numerical 
accuracy, in identical energies and densities. Furthermore, for vanishing barrier height, i .e., for A = 0, our solutions 
for cylindrically symmetric harmonic traps agree with those reported in the literature [la . l54j . For non- vanishing 
barrier height, we compared our solutions with those reported in Ref. [!§]. Our energies and chemical potentials are 
in reasonable agreement with those reported in Ref. (49|. For A — 4E Z , b — 0.2a z , A = 0.1 and D = 0.6, e.g., we 
find E/N — 10.69£ z and = 10.00i? z while the values reported in Fig. 2 of Ref. f4§] are smaller by about 3 and 4%, 
respectively. These deviations are somewhat larger than our estimated numerical uncertainty. 

The time dynamics of the system is determined by evolving a given initial state in real time. The initial state 
is chosen according to the variational two-mode model wave function (see Sec. Illlj) or by adding a small random or 
smooth perturbation to the stationary GP wave function of the energetically lowest lying state. If the system collapses 
to a high density state in response to the application of a small perturbation, then our simulations are only able to 
follow the real time evolution for a limited time period. Eventually, our grid becomes too coarse to accurately present 
the time evolved state. Since our main aim is directed at identifying the stability and the collapse mechanisms, this 
artefact does not pose any true limitations on our analysis. In fact, once the density becomes sufficiently high, the 
mean-field GP description breaks down anyways and beyond mean-field corrections need to be included. Such a 
treatment is, however, beyond the scope of the present work. 



B. Bogoliubov de Gennes equations 



In addition to time-evolving a given initial state, we analyze the stability of the dipolar BEC by seeking solutions 
to the time-dependent GP equation of the form [56[ 

1>{f, t) = cxpi-ifit/h) [V-o(r) + 5i/>{f, t)] , (11) 

where tpo(r) denotes the energetically lowest lying solution of the time- independent GP equation with k = and /i 
the corresponding chemical potential. We seek "perturbations" Sip(r,t) that oscillate with frequency ui, 

5ip(r,t) = u(r) exp(— iuit) + v* (r) exp(iurf), (12) 

where u(r) and v(r) denote the Bogoliubov de Gennes "particle" and "hole" functions Plugging Eq. (jTTJ) with 
Sip given by Eq. (|12p into Eq. fTJ), keeping terms up to first order in 6i/j(r, t) and its complex conjugate, and equating 
the coefficients of the terms oscillating with exp(— iut) and exp(iwi), respectively, we find the Bogoliubov de Gennes 
equations [5J] 

hum(r) — A(r)u(r) + 
(N-l)J V/ dd (r~r^oVM^)dV^ (r) + 

{N-l)[ V dd (r-^)Mr)v(r)d 3 ^Mr) (13) 

and 

— Hujv*(f f ) — A(f)v*(r) + 
(N-l)J V AA {r-?)r [?)v* Mr) + 
(N - 1) J V dd (?- ^)V>o(^V)tfWo(rO- (14) 
In Eqs. (fT3|) and (fl"4"]) . the operator A(r) is defined as 

A{r) = H -V + 

(N - 1) / V dd (r- ^I^^IW, (15) 



where Hq denotes the Hamiltonian of the non-interacting system, 

n 2 



H a = -^—V 2 + V CKt (r). (16) 
Ira 
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Equations (fT3f and (fT4|) can be decoupled by introducing two new functions / and g, f(r) 
j(f) = —u{r) + v(r). Assuming, without loss of generality, that ipo(r) is real, we find 



2{N-l)A{r) 



h 2 u; 2 f(r) = A(f) [A(f)f{r)\ + 
/(r')Vdd(f - 7)M7)d^Mr 



u(r) + v(r) and 



(17) 



and 



h 2 u 2 g{r) = A{r) [A(f)g(r)} + 
2(N -I) J Vd d (f - ^)M^)M^)9^)d^M^- (18) 

Following Ref. [54| . we solve Eq. (fl8|) for the square of the Bogoliubov de Gennes excitation frequency ui and the 
corresponding eigenvector g(f) iteratively using the Arnoldi method. Once g(r) is determined, the eigenvector f(r) 
can be obtained from the identity 



m = -^A(f)g(f). 



(19) 



The physical meaning of / is elucidated by calculating the density \ij>(r, t)\ 2 up to first order in Sip and its complex 
conjugate. For real u and v, this gives 



(20) 



which shows that /(f), together with ipo(^) and u>, determines the time-dependent density. Due to the cylindrical 
symmetry of the system, the ip dependence of /(f) separates, /(f) = f(p, z)h((p). Section [V] discusses the behavior of 
/(/?, z), which we refer to as the Bogoliubov de Gennes eigenmode, for different k and various (D, A) combinations. 

The outlined approach allows for the determination of a sequence of excitation frequencies for a given azimuthal 
quantum number k at a time. It can be seen from Eq. (|12p that a negative u> 2 and thus a purely imaginary ui 
corresponds to a situation where the stationary ground state solution is dynamically unstable. 



III. TWO-MODE MODEL 



Atomic BECs, coupled through non-vanishing potential barriers, have been used extensively to model coupled 
condensed matter systems such as 3 He-B reservoirs [HI, Hi], |42|, |45|, [58[. Although neutral, the study of weakly- 
coupled atomic BECs allows, e.g., for the realization of a variety of typical dc and ac effects that characterize charged 
Cooper pair superconducting junctions [42l . [59| . The connection between weakly-coupled atomic BECs and more 
traditional condensed matter systems becomes most apparent if the former is approximated by a two-mode model 
and mapped to a Josephson like Hamiltonian. Here, our primary motivation for employing the two-mode model is to 
develop an intuitive understanding of some of the phenomena observed in our time-independent and time-dependent 
mean-field studies. 

Let ipsir) and ipA(r) denote the energetically lowest lying stationary GP solutions that are respectively symmetric 
and anti-symmetric with respect to z = 0. If the symmetric function ipsi'r) is the energetically lowest lying solution of 
the stationary GP equation, we calculate it by employing the conjugate gradient method or by evolving in imaginay 
time (see Sec. |TTJ) . The anti-symmetric solution ^a(0 is obtained by restricting the basis functions employed in the 
conjugate gradient method to functions that are anti-symmetric with respect to z = 0. Without loss of generality, we 
assume in the following that ips and ipA are real. In the two-mode model, the solutions ifis and i/'A are treated as a 
basis that defines the two "modes" $l(0 and <!>R,(r), 

* £ .*(f> = ^ S(r1± ^ A(f0 - ( 21 ) 
V2 

By construction, $i(f) and $i?(f) are normalized to one and orthogonal to each other. The functions $l and <1>r 
are, for appropriately chosen parameters, located predominantly in the left well and in the right well, respectively. 
Within the two-mode model, the time-dependent wave function is approximated by (see, e.g., Ref. [5.(5:]) 



ip(r,t) = c L (t)$ L (r) + c R (t)$ R (r), 



(22) 
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where the complex- valued time-dependent expansion coefficients cl(£) and c R (t) are related through the normalization 
condition |cl(£)| 2 + I c r(*)| 2 = 1- Defining c L . R {t) — \cl,r{P}\ ex p[*^L,_R(i)], the time evolution within the two-mode 
model is governed by two variables: the fractional difference Z(t) of the population located in the left and in the right 
well, 

Z(t) = \c L (t)f-\c R (t)f, (23) 
and the phase difference or relative phase <f)(t), 

#t) = 0a(*) -&-(*)• ( 24 ) 

Plugging Eq. (f2"2"|) into Eq. ([I}, multiplying by and <I> R (r), respectively, and integrating out the spatial degrees 

of freedom, we obtain two coupled equations that govern the time dynamics: 

. fc <M*) 

in 

fir 

[E + B + (U- B)\c L {t)\ 2 ] c L (t) - Tc R (t) (25) 

and 

. fc dc R (i) 

in : — — 

dt 

[E + B + (U- B)\c R {t)\ 2 ] c R (t) - Tc L (t). (26) 
In deriving Eqs. ([25]) and (f26|) . we neglegted terms of the form 

(N — 1) X 

M^j(r)V dd (r- (27) 



with i ^ j or k ^ I, where k and I can take the values L and R. These terms are small as long as $l and <I> R are 
located predominantly in the left well and in the right well, respectively. In Eqs. (|25|) and (|26|) . the onsite, offsite (or 
interaction tunneling) and tunneling matrix elements U, B and T are defined as 

U=(N-l)x 

($ L (r1)V dd (f- fO($ L (rO)Wd 3 f, (28) 



B = (N- 1) x 



/ J ($ L (f)) 2 F dd (f'-r')($ R (r')) 2 ( i 3 r'd 3 f 



(29) 



and 



T 



— V$ L (r)-V$ R (rO- 
2m 

$ L (r)Kxt(r)* R (^]d 3 r, 



and the "zero point energy" i^o is defined as 



En — 



— |V$ L (r1| 2 + y oxt (r)(<i> L (r)) 2 
2m 



(30) 



(31) 



Usage of 3> R (r) instead of Q^if) in Eqs. (|28|) and (l3T|) gives the same result. 

Rewriting the coupled equations (|23|) and (|26p in terms of Z(t) and 0(t) leads to the classical Hamiltonian Htm 
(using h = 1), 



(32) 
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FIG. 1: (Color online) Character of the energetically lowest lying stationary CP solution for b — 0.2a z and A — 12E Z : The 
dash-dot-dotted, dash-dotted and dashed lines indicate those A and D values at which the character of the energetically lowest 
lying stationary CP solution with k — changes from symmetric (S) to symmetry-broken (SB), from symmetry-broken to 
unstable (U), and from symmetric to unstable, respectively. The "solid (red) islands" in the upper right corner of the phase 
diagram (comparatively large D and A) indicate two regions of the phase diagram where the solutions are symmetric but where 
the density maximum is located at p > 0; these islands are discussed in more detail in the context of Fig. [5] For comparison, 
circles show the boundary between the symmetric and symmetry-broken regions predicted by the two mode model. Note the 
log scale of both axes. 



where 

A=^. (33) 

Notably, Z(t) and 4>(t) are conjugate variables of the classical Hamiltonian. The energy of Htm is conserved and 
can, e.g., be obtained by inserting Z(0) and </>(0) into Eq. (f3"2")l . The properties of Htm have been discussed in detail 
in the literature [39L l4ll l42l. |45| . Here, we review a few points that will aid in the understanding of our GP solutions. 

The two-mode model immediately leads to three different classes of stationary solutions, i.e., solutions with constant 
Z(t) and 4>(t): A symmetric solution for <j>(t) = 2im (n integer) and Z(t) = 0; its energy is —2T. An anti-symmetric 
solution for <fr(t) = (2n+l)7r (n integer) and Z(t) = 0; its energy is 2T. A symmetry-broken solution for (j>(t) = (2n+l)-7r 
(n integer) and Z(t) — ±vl — A -2 ; this solution exists only if |A| > 1 and its energy is T(A + A -1 ). Section HVl 
compares these stationary two-mode model solutions with those obtained from the stationary GP solutions. 

It turns out that Hamilton's equations of motion can be solved analytically for Htm [HI]. Of particular interest for 
our study is the so-called Josephson oscillation frequency luj, which — for small amplitude motion — can be expressed 
in terms of A (with h "restored" ) , 

2TVTTX (34) 

Section IVl compares the two-mode model frequency cjj^tm with the frequency obtained from the real time dynamics 
and by solving the Bogoliubov de Gennes equations. For the real time dynamics, we prepare an initial state at time 
t = according to Eq. (|2"2"|) and then time evolve this state according to the time dependent mean-field Hamiltonian. 
A Fourier analysis of the expectation value of z(t) then reveals the predominant excitation frequency. 



IV. DISCUSSION OF STATIONARY SOLUTIONS 



This section discusses our solutions to the stationary GP equation. In particular, we present the phase diagram as 
a function of the aspect ratio A and the mean-field strength D for a fixed barrier height A and discuss selected density 
profiles. Furthermore, we discuss how the phase diagram changes with varying barrier height and explain some of the 
GP results within the two-mode model. 

Figure Q] summarizes the character of the energetically lowest lying solutions with k = of the stationary GP 
equation as functions of the aspect ratio A and the mean- field strength D for A = 12E Z . The parameter combination 
(A, D) = (0.1,1) corresponds, e.g., to a Cr condensate with vanishing s-wave scattering length, oj z — 2n x 10Hz, 
ojp = 2tt x 100Hz and N w 1835. The "phase diagram" consists of three regions: Firstly, a region where the 
energetically lowest lying state with k = of the stationary GP equation is symmetric with respect to the z-axis; we 
refer to this solution as symmetric ("S") throughout this paper. Examplary density profiles are shown in Figs.J^a), 
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FIG. 2: (Color online) Density plots of the energetically lowest lying stationary GP solution with k = for b = 0.2a z , A — 12E Z 
and four different (D,\) combinations: (a) {D,\) = (0.5477,0.3) (symmetric solution), (b) {D,X) = (1.643,0.3) (symmetry- 
broken solution), (c) (D, A) = (316.2, 10) (symmetric solution with density maximum at p = 0) and (d) (D,X) — (474.3, 10) 
(symmetric solution with density maximum at p > 0). The contour lines are chosen equidistant in all four panels. The dashed 
contours correspond to (a) 0.05aJ 3 , (b) O.laJ 3 , (c) O.OOOlaJ 3 and (d) O.OOOlaJ 3 , while the solid contours correspond to (a) 
0.35aJ 3 , (b) 0.7aJ 3 , (c) 0.0007aJ 3 and (d) 0.0007aJ 3 . 



[He) and HJd) (see below for more details). Secondly, a region where the energetically lowest lying state of the 
stationary GP equation is neither symmetric nor anti-symmetric with respect to the z-axis; we refer to this solution 
as symmetry-broken ("SB") or asymmetric. An examplary density profile is shown in Fig. HJb). And thirdly, a region 
where the stationary GP equation supports a high-density or collapsed solution but no gas-like solution. We refer to 
this solution as mechanically unstable ("U"). Sections IIV Al and IIV Bl discuss the properties of the phase diagram in 
more detail for A < 1 and A > 1, respectively. 



A. "Small" aspect ratio (A < 1) 

Figure [21a) shows the energy contributions to the total energy per particle E tot /N for A = 0.3 and A = 12E Z as a 
function of D: the kinetic energy per particle i?kin (dashed line); the harmonic trap energy per particle i?ho (dotted 
line), which is defined as the expectation value of Vh \ the Gaussian energy per particle E g (dash-dotted line), which 
is denned as the expectation value of V g ; and the mean- field dipole-dipole energy per particle E'dip (dash-dash-dotted 
line), which is defined as the expectation value of the mean-field term [third term on the right hand side of Eq. ©]. 
A solid line shows the total energy per particle E tot /N. The energies terminate at the critical value D CI at which the 
stationary GP equation first supports a negative energy solution. 
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FIG. 3: (Color online) (a) Energy contributions of the energetically lowest lying stationary GP solution with k = as a function 
of D for A = 0.3, A = 12E Z and b — 0.2a z . The solid, dashed, dotted, dash-dotted and dash-dash-dotted lines show Etot/N, 
£-kin, £"ho, E g and -Edip, respectively, (b) Blow-up of the Gaussian energy E s . E s exhibits a kink at D ~ 0.75, indicating the 
symmetry change (symmetric to symmetry-broken) of the energetically lowest lying stationary GP solution. 



The Gaussian energy E g is shown on an enlarged scale in Fig. OJb). It can be seen that E g shows a "kink" at 
D w 0.75. We find that the other energy contributions (i.e., E]& n , Eho and -Edip) and E tot /N exhibit kinks at the same 
D value. These kinks are, however, less pronounced and not (or hardly) visible on the scale shown in Fig. |3fa). Our 
analysis shows that the D values at which the kinks occur coincide with the D values at which the density profiles 
of the energetically lowest lying stationary GP solutions change from symmetric to symmetry-broken. In most of 
our calculations for fixed 6, A and A but varying D, we use the kink in E g to determine the D value at which the 
character of the energetically lowest lying stationary GP solution changes and thus to obtain the dash-dot-dotted line 
in Fig.[TJ The stability of the solutions around the symmetry to symmetry-broken transition is discussed in Sec. [V] in 
the context of Figs. 171 through fTUl 

The dash-dot-dotted lines in Fig. Q] can be reproduced qualitatively by the two-mode model (see circles in Fig. [J). 
To illustrate some aspects of the two-mode model, solid and dashed lines in Fig. 2] show A [see Eq. ([55]) ] as a function 
of D for A = 0.3 and 0.4, respectively, and A — 12E Z and b = 0.2a z . For |A| < 1 and positive T, the two-mode model 
predicts a symmetric stationary ground state. For |A| > 1, a symmetry-broken solution is supported; if T > and 
A < — 1, the symmetry-broken state has a lower energy than the symmetric state. Vertical arrows in Fig. [4] mark the 
D values, D w 1.31 and 2.89, at which the transition from symmetric to symmetry-broken occurs for A = 0.3 and 
A = 0.4, respectively. These two-mode model predictions (also shown as circles in Fig. [I} are slightly larger than the 
results obtained by solving the GP equation but predict the symmetric to symmetry-broken transition qualitatively 
correctly. 

It is interesting to compare the behavior of A, which can be interpreted as the ratio between an effective interaction 
energy and twice the tunneling energy, for A = 0.3 and 0.4 (solid and dashed lines in Fig. [|]). For A = 0.3, an increase 
of the mean-field strength D leads to a monotonic decrease of A. For A = 0.4, in contrast, A first increases, reaches 
a maximum at D ss 0.87 and then decreases monotonically. We find that the onsite energy U and the offsite energy 
B are both negative for all D shown in Fig. [U A change of the aspect ratio A effectively changes the strength of the 
dipolc-dipole interaction, leading to a less attractive U than B, and thus to a positive A, for small D and A = 0.4. For 
A = 0.3, in contrast, the onsite energy U is always more negative than the offsite energy B, resulting in a negative A 
for all D. 

In addition to the barrier height A — Y1E Z , we considered smaller barrier heights A, in particular A = AE Z and 8E Z , 
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FIG. 4: (Color online) Two-mode model parameter A for A — 12E Z and b = 0.2a z as a function of D for A = 0.3 (solid line) 
and A = 0.4 (dashed line). Vertical arrows mark the D values at which |A| equals 1; for |A| < 1 and > 1, the two-mode model 
predicts that the energetically lowest lying stationary state is symmetric and symmetry-broken, respectively. 

for a few selected A values. Our calculations suggest that the dash-dot-dotted line in Fig.[T](i.e., the line that marks the 
symmetric to symmetry-broken transition) moves to larger D values with decreasing A while the dash-dotted line (i.e., 
the line that marks the symmetry-broken to unstable transition) remains approximately unchanged with decreasing A. 
The dependence of the dash-dot-dotted line on A for fixed A and b can be explained by applying the two-mode model. 
As A decreases, the tunneling energy T becomes more important compared to the absolute value of the effective 
interaction energy U — B. This implies that |A| decreases with decreasing A (for fixed A and b). Correspondingly, 
a larger D is required for the two-mode model condition |A| = 1, which signals the symmetric to symmetry-broken 
transition, to be fulfilled. The fact that the dash-dotted line in Fig. Q] remains to first order unchanged with decreasing 
A is due to the fact that the density of the system prior to collapse is located predominantly in one of the wells. This 
implies that the density prior to collapse is only weakly dependent on A, thus explaining the comparatively small 
dependence of the dash-dotted line on A for the parameter combinations investigated. 

We note at this point that the linear stationary Schrodinger equation permits only symmetric and anti-symmetric 
solutions but no symmetry-broken solutions. This fact emphasizes that the transition from symmetric to symmetry- 
broken is driven by mean-field interactions. Furthermore, this fact implies that the symmetry-broken solution should 
disappear if sufficiently many higher order corrections to the mean-field GP equation are taken into account (see, 
e.g., Ref. [41]). In this sense, the appearance of the symmetry-broken region in the phase diagram is an artefact of 
the mean-field formalism. It is, however, intimately related to the dynamical phenomena of Josephson oscillation and 
macroscopic quantum self-trapping, both of which have been observed experimentally for s-wave interacting BECs. 
We return to these considerations in Sec. |V| m the context of the discussion of Figs. [71 through [TU1 

B. "Large" aspect ratio (A > 1) 

Figure [5] shows an enlargement of the large A region of Fig. Q] using a linear scale for both A and D. The So region of 
the phase diagram is characterized by GP solutions whose density maxima are located at p — [see Figs.[5Ja) and[^c) 
for examples] while the S>o region of the phase diagram is charactericed by GP solutions whose density maxima are 
located at p > [see Fig.[2][d) for an example]. The latter class of density profiles only exists in a narrow parameter 
region of the phase diagram; in particular, these solutions only arise for pancake-shaped confining potentials and not 
for cigar-shaped confining potentials. Furthermore, the solutions with S>o character are unique to dipolar gases, i.e., 
they are not observed for purely s-wave interacting gases, and thus directly reflect the anisostropic long-range nature 
of the dipolc-dipole interactions. 

The two different classes of symmetric solutions have previously been characterized for A = 0, i.e., for a pancake- 
shaped trapping geometry without barrier [TH, |60T]. In those studies, a dipolar BEC with density maximum at p > 
was termed "red blood cell" , as its isodensity surface is reminiscent of the shape of a red blood cell. The S>o regions 
in Fig. are characterized by the formation of two staggered red blood cells. Section [V] shows that the dynamical 
instability near D cr of the stationary k = ground state solutions of types So and S>o is distinctly different. 

Figure [5] shows that, generally speaking, the D value at which the dipolar gas becomes unstable increases with 
increasing A. This trend can be understood by realizing that an increase of A leads to a "flattening" of the system so 
that the dipoles interact effectively more repulsively. The boundary near the stable and unstable regions shows a rich 
structure: i) As already noted above, S>o islands in which the density profiles are structured exist, ii) The boundary 
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FIG. 5: (Color online) Blow-up (with more detail) of the phase diagram for b — 0.2a z and A = 12E Z shown in Fig. [T] The 
region where the energetically lowest lying symmetric stationary GP solution has its density maximum at p — is labeled by 
"So" and that where the energetically lowest lying symmetric stationary GP solution has its density maximum at p > by 

B>0 ■ 



between the So and the U regions of the phase diagram changes non-monotonically. For D w 240, e.g., the system is 
mechanically stable for A > 8.13, mechanically unstable for 8.13 > A > 7.72, and then again mechanically stable for 
a small A regime (7.72 > A > 7.42). 

We find that some, though not all, of the features of the phase diagram can be reproduced qualitatively by a simple 
variational wave function ip va .r(p, z), 
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where b\ — bg denote variational parameters that are optimized by minimizing the energy per particle. For 62 = 65 = 0, 
ip V ar reduces to the commonly used variational wave function of purely Gaussian shape. The second term in the first 
square bracket on the right hand side of Eq. (|35l) has been added to allow for the description of densities of S>o 
character while the second term in the second square bracket on the right hand side of Eq. (|33|) has been added to 
account for the Gaussian barrier along the z-direction. Figure [5] compares the total energy per particle from our 
variational calculation (dashed line) with that from the full numerical calculation (solid line) for A = 12E Z , A = 7 and 
b = 0.2a z . The variational energy is less than 2 % higher than the energy obtained from the full numerical calculation. 
We find that the density of the dipolar gas changes from So to S>o character at D ss 35, compared to D = 80.03 
obtained from the full calculation. For both sets of calculations, the energy and its derivative change smoothly as 
the system undergoes the structural change from So to S>o- For comparison, a dash-dotted line in Fig. [6] shows the 
energy per particle for ip va r with 62 = (we refer to this variational wave function as four-parameter wave function) , 
i.e., for a wave function that is not sufficiently flexible to describe structured ground state densities of red blood cell 
shape. As expected, this variational wave function results in somewhat higher energies. 

The variational wave function ip var predicts So to S>o transitions for all aspect ratios A between 5 and 12, indicating 
that it is not flexible enough to describe the island character of the S>o regions of the phase diagram and, furthermore, 
that the So to S>o transition is driven by a delicate balance between the different energy contributions. Motivated 
by calculations presented in Ref. (l5j . we expect that the variational four-parameter wave function can qualitatively 
reproduce the existence of alternating stable and unstable regions of the phase diagram as A is changed for fixed D 
and A (see our discussion above for D w 240 and A = 12E Z ); we have, however, not checked this explicitly. Lastly, 
we note that the variational six-parameter wave function predicts a stable dipolar gas even for fairly large D (i.e., D 
values larger than those shown in Fig. [6]) while the full numerical calculation predicts collapse at D w 190.49. 



V. DISCUSSION OF DYNAMICAL STUDIES 



This section presents Bogoliubov de Gennes excitation spectra and discusses the corresponding eigenmodes. For 
small A and appropriate D (see Sec. IV A[) . the lowest non- vanishing Bogoliubov de Gennes excitation frequency is 
identified as the Josephson oscillation frequency luj. Comparisons with results obtained by time-evolving a properly 
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FIG. 6: (Color online) A solid line shows the total energy per particle Etot/N of the energetically lowest lying stationary GP 
solution with k = as a function of D for A = 7, A — 12E Z and b = 0.2a z obtained numerically. For comparison, dashed and 
dash-dotted lines show Etot/N obtained using the variational six- and four-parameter variational wave functions (see text for 
details). The density of the system changes from So to S>o character at D w 35 and 80.03 for the six-parameter variational 
and the full numerical calculations, respectively. 




FIG. 7: (Color online) Excitation spectrum obtained by solving the Bogoliubov de Gennes equations as a function of D for 
A = 12E Z , b = 0.2a z and A = 0.3. The real parts of the frequencies for k = and 1 are shown by solid and dashed lines, 
respectively. The vertical arrow indicates the D value, D ~ 0.68, at which the real part of the lowest non-vanishing k = 
Bogoliubov de Gennes frequency vanishes. 



prepared initial state and by applying the two-mode model equations are presented. In the regime where the sym- 
metric solutions are of types So and S>o, respectively (large A, see Sec. IV Bp . the decay mechanisms are identified. 
Furthermore, the character of various (avoided) crossings of the excitation frequencies is revealed. 



A. "Small" aspect ratio (A < 1) 

Figure [7] shows the excitation spectrum as a function of D obtained by solving the Bogoliubov de Gennes equations 
for A = 12E Z , b — 0.2a z and A = 0.3. The spectrum is characterized by three distinct features that will be elaborated 
on in the following paragraphs: i) The real part of the lowest k = frequency vanishes at D ss 0.68, and "reappears" 
at D w 0.75. ii) The k = frequencies show a series of crossings (or avoided crossings) at D sa 2.42. iii) At slightly 
larger D values, i.e., near D w 2.45, the real part of several k = excitation frequencies vanishes. 

We first discuss the regime i) around D ps 0.68—0.75. Figure[5]Ja) shows the Bogoliubov de Gennes eigenmode f(p, z) 
that corresponds to the lowest non-vanishing k = frequency for D = 0.6573, A = 12E Z , A = 0.3 and b — 0.2a z . 
For these parameters, the energetically lowest lying stationary GP solution is symmetric and the corresponding 
eigenfrequency has a finite real part and vanishing imaginary part (see Fig. [7]). Since Bogoliubov de Gennes functions 
with k — have no explicit tp dependence, the eigenmode shown in Fig. [SJa) corresponds to a situation where the 
population oscillates with frequency to between the left and the right well as a function of time. The lowest non- 
vanishing k = frequency can thus be identified as the Josephson oscillation frequency toj (see also below). For 
comparison, Fig. [SJ^b) shows the Bogoliubov de Gennes eigenmode f(p,z) corresponding to the lowest non- vanishing 
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FIG. 8: (Color online) Bogoliubov de Gennes eigenmodes f(p,z) corresponding to the lowest non-vanishing k = frequency 
for A = 0.3, A = 12E Z , b = 0.2a z and (a) D = 0.6573 and (b) D = 0.7668. The contours are chosen equidistant, with solid and 
dashed lines corresponding to positive and negative values of /. The dash-dotted lines indicate the nodal lines of /. 



k = frequency for D — 0.7668 (i.e., in the regime where the frequency has "reappeared" and where the energetically 
lowest lying stationary GP solution with k = is symmetry-broken) and the same A, E z and b values as before. In 
this case, the asymmetry of the eigenmode indicates that there is population transfer between the left and the right 
wells but that there is, on average, more population in the right than in the left well. This behavior is identified as 
macroscopic quantum self-trapping. Our interpretation of the Bogoliubov de Genne eigenmodes is supported by our 
time-dependent calculations. 

In our time-dependent studies near the symmetric to symmetry-broken transition, we prepare an initial state and 
time evolve it according to the mean-field Hamiltonian H, Eq. As for s-wave interacting BECs, the system 

dynamics can be divided into two categories (see also above): A regime where the population is transferred back and 
forth between the left well and the right well (this is the Josephson oscillation regime) and a regime where the time 
averaged population is asymmetrically divided among the two wells (this is the macroscopic quantum self-trapping 
regime). Figures[^a) and[5Jb) show the time evolution of the expectation value (z(t)) for D = 0.6573 and D — 0.7668, 
respectively. Here, (z(t)) is obtained by calculating the expectation value of z with respect to the GP density at each 
time step. The expectation value (z(t)} is related to but not identical to the population difference Z(t) introduced in 
Sec. IIIII For these D values, the energetically lowest lying stationary GP solution is symmetric and symmetry-broken, 
respectively. For D = 0.6573, the initial state is prepared according to Eq. with (f)(0) = and Z(0) = 0.002. 
Figure Ufa) shows that (z(t)) oscillates between positive and negative values of equal magnitude and that the time 
average of (z(t)) over a period gives zero. For D = 0.7668, the initial state is prepared by adding a small amount of 
random noise to the energetically lowest lying stationary GP solution. Figure [S^b) shows that (z(t)) oscillates about 
a negative value and that the time average of (z(t)} gives a non-zero value. An analysis of the time evolution of the 
density profiles confirms that the system is in the Josephson regime and in the macroscopic quantum self-trapping 
regime, respectively. 

To determine the oscillation frequency from the time evolution of (z(t)), we Fourier transform (z(t)) and record 
the center of the dominant peak for various parameter combinations. Circles in Fig. 1101 show the resulting Josephson 
oscillation frequency Wj for A = 12E Z , X = 0.3 and b = 0.2a 2 . The agreement between the frequency obtained from the 
Fourier analysis (circles in Fig. [TTJ)l and the lowest non-vanishing k = Bogoliubov de Gennes excitation frequency 
(solid line in Fig. [TTJ)l is excellent. The D value at which the Josephson oscillation frequency obtained by Fourier 
transforming (z(t)) vanishes, coincides, within our numerical accuracy, with that at which the lowest non- vanishing 
k = Bogoliubov de Gennes excitation frequency becomes imaginary. Notably, this D value, D 0.68, is slightly 
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FIG. 9: (Color online) Solid lines show the expectation value (z(t)) (see text) — calculated by time evolving a given initial state 
according to the mean-field Hamiltonian, Eq. ((2} — as a function of time t for A = 12E Z , A = 0.3, b = 0.2a z and (a) D — 0.6573 
(Josephson oscillation regime) and (b) D = 0.7668 (macroscopic quantum self-trapping regime). In panel (a), the initial state 
is prepared according to Eq. (|22p with Z(0) = 0.002 and 0(0) = 0. In panel (b), the initial state is prepared by adding a small 
amount of random noise to the energetically lowest lying stationary GP solution. 




FIG. 10: (Color online) Josephson oscillation frquency ujj as a function of D for A = 0.3, A — 12E Z and b — 0.2a z . The circles 
show the Josephson oscillation frequency uj obtained from our time-dependent study, in which the initial state is prepared 
according to Eq. (|22|1 with Z(0) = 0.002 and 0(0) = and then time evolved according to the mean-field Hamiltonian H , 
Eq. The solid line shows the lowest non- vanishing k — Bogoliubov de Gennes excitation frequency. For comparison, a 
dash-dotted line shows the two-mode model prediction luj.tm- 



smaller than the D value at which the energetically lowest lying stationary GP solution changes from symmetric to 
symmetry-broken [D « 0.75). 

We find that the lowest non-vanishing k = Bogoliubov de Gennes frequency for D — 0.7668 and A = 0.3 is about 
30% smaller than the oscillation frequency extracted from Fig. Efb), i.e., oj = 0.059ijJ 2 . The fact that the Bogoliubov 
de Gennes excitation frequency differs notably from the frequency obtained by Fourier-transforming (z(t)) might be 
due to the approximate nature of the Bogoliubov de Gennes equations. 

For comparison, a dash-dotted line in Fig. [TTJ1 shows the Josephson oscillation frequency u> j,tm, Eq. (|34|>. predicted 



15 




c3 



C3 



0.4 0.8 1.2 



3 

-3 





(c) : 










. i . . 





0.4 0.8 1.2 



FIG. 11: (Color online) Bogoliubov de Gennes eigenmodes f(p,z) corresponding to the four lowest non-vanishing k — 
frequencies for A = 0.3, A = 12E Z , b = 0.2a z . Panels (a)-(d) show / corresponding to the lowest, second lowest, third lowest 
and fourth lowest frequencies for D — 2.410 (i.e., just before the crossing) while panels (e)-(h) show / corresponding to the 
lowest, second lowest, third lowest and fourth lowest frequencies for D = 2.443 (i.e., just after the crossing). The contours are 
chosen equidistant, with solid and dashed lines corresponding respectively to positive and negative values of /. 



by the two-mode model. Figure [TO] shows that the two-mode model provides a qualitatively but not quantitatively 
correct description of the Josephson oscillation frequency. The fact that the two-mode model does not allow for 
quantitative predictions for all D is likely due to the fact that the modes $l and <I>r, are not entirely located in the 
left well and in the right well, respectively, but that the left mode "leaks" into the right well and the right mode into 
the left well. This has been discussed in some detail in Ref. [4!|, which employs a slightly modified version of the 
two-mode model. In an attempt to obtain a better simple quantitative description of the system dynamics, we applied 
the improved two- mode model proposed in Ref. (45|. For the cases considered, we find that this model leads only to 
small changes compared to the simple two-mode model applied above and does not provide a significantly improved 
description. In the future, it may be interesting to apply a multi-mode model. 

We now discuss the regime ii) near D « 2.42, where the k = frequencies show (avoided) crossings. To shed light 
on these (avoided) crossings, Figs. lllf a)-(d) show the Bogoliubov de Gennes eigenmodes f(p, z) corresponding to the 
four lowest non-vanishing k = frequencies just before the crossing (i.e., for D = 2.410), while Figs. HTT e)-(h) show 
those corresponding to the four lowest k = frequencies just after the crossing (i.e., for D = 2.443). Dash-dotted 
lines in Fig. QT] indicate the nodal lines of the Bogoliubov de Gennes eigenmodes. While some of these nodal lines are 
to first order only dependent on z, others depend in a non-trivial manner on p and z. In the following we discuss a 
few key features of the eigenmodes shown in Fig. 1111 The eigenmode corresponding to the lowest frequency extends 
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over both wells just before the crossing [see Fig. flTT a)] and is located predominantly in one of the wells just after 
the crossing [see Fig. Illf e)]. The eigenmode corresponding to the second lowest non- vanishing frequency, in turn, is 
located predominantly in one of the wells just before the crossing [see Fig. [TTTb)] and extends over both wells just 
after the crossing [see Fig. flTTf)]. The third and fourth lowest excitation frequencies show an avoided crossing at 
D « 2.41 (see Fig. [7]). The eigenmode corresponding to the third lowest non- vanishing k = frequency has a fairly 
small amplitude in the right well before the avoided crossing [see Fig. fTTl c)]: after the avoided crossing, the eigenmode 
corresponding to the third lowest k = frequency is located essentially entirely in the left well [see Fig. ITIf g)]. The 
eigenmode corresponding to the forth lowest frequency, in turn, changes comparatively little as D increases [see 
Figs. \H\d) and filth)]. 

Next, we turn to regime iii) near D 2.45, where the real part of several k — Bogoliubov de Gennes ex- 
citation frequencies vanishes. The softening of these frequencies is inherently related to the transition from the 
stable symmetry-broken region to the unstable region of the phase diagram. The Bogoliubov excitation spectrum for 
A z = 12E Z , A = 0.3 and b — 0.2a z shows that the real part of the lowest k = frequency vanishes at D ss 2.45, which 
is slightly smaller than the D value at which the stationary GP equation starts supporting unbounded negative energy 
solutions. In particular, for the parameter combination considered, the difference is about 0.6%. The Bogoliubov de 
Gennes excitation spectrum indicates that the collapse is triggered by the lowest k = mode. The corresponding 
eigenmode [see Fig. Illf e)] shows that, as might be expected naively, the density grows appreciably in the well that 
supports the majority of the population. This interpretation is supported by our time-dependent studies. Following 
the lowest k = mode, the eigenmode corresponding to the third lowest non-vanishing k = frequency becomes 
soft. As indicated in Fig. [TlTg), this eigenmode is, not unexpectedly, also located predominately in one of the wells. 
To summarize, the collapse of the system can be characterized as a global collapse in which the density maximum 
increases in one of the wells, with the density maximum being located at p = 0. 

B. "Large" aspect ratio (A > 1) 

This section discusses selected Bogoliubov de Gennes excitation spectra for larger A. In particular, we focus on 
A = 5 and A — 7, for which the energetically lowest lying stationary GP solution prior to collapse is of type So and 
5>o, respectively. 

Figure [12] shows the excitation spectrum as a function of D for A = 5, A = 12E Z and b = 0.2a z . The two lowest 
non-vanishing k — frequencies [solid lines in Fig. [T2Ta)] cross at D w 17. Prior to the crossing, the lowest k = 
frequency corresponds to an eigenmode whose nodal line is parameterized by z = [see dash-dotted line in Fig. ll3f a)]. 
Correspondingly, the eigen mode describes density oscillations between the left and the right well with, on average, 
equal densities in each of the two wells. For D 17—22 (i.e., after the crossing), the nodal line of the lowest 
non-vanishing k = eigenmode is to a good approximation independent of z and can be parametrized by p « 2a z 
to « 3a z . As an example, Fig. [T37 b) shows the eigenmode for D = 21.87, which is just slightly smaller than the 
D value at which the real part of the corresponding Bogoliubov de Gennes frequency vanishes. In particular, the 
lowest Bogoliubov de Gennes mode becomes soft at D w 21.92 while the stationary GP equation ceases to a support 
a positive energy solution with k — at a somewhat larger D value, namely at D w 22.03. 

Figures [12] and [TBT b) suggest that the collapse of the system is triggered by the lowest k — mode and that the 
collapse is associated with an increase of the peak density at p = in each of the two wells. The collapse can thus 
be characterized as a local collapse as opposed to a global collapse. The local nature of the collapse (i.e., the fact 
that the peak density grows simultaneously in two distinct regions) can be traced back directly to the presence of 
the comparatively large Gaussian barrier. The chemical potential p takes values around 2E Z for the parameter range 
considered in Fig. [T2l and is thus significantly smaller than the barrier height A, A = 12E Z . In the limit that the 
Gaussian barrier vanishes [ljl, the collapse becomes global. In this case, a similar nodal pattern of the eigenmode 
was found (see Fig. 21c of Ref. [l5[ ; the larger number of nodal lines in this plot can be traced back to the larger A 
value) and the collapse is associated with a radial roton. 

Figure [T4l shows the Bogoliubov de Gennes excitation spectrum for A = 7, A = 12E Z and b = 0.2a z . For this 
parameter combination, the energetically lowest lying stationary GP solution deviates from the "structureless Gaussian 
shape" prior to collapse and is instead characterized by a density whose maximum is located at p > 0, i.e., the density 
is of type S>o prior to collapse (see Fig. [2] for a density profile of type S>o for a somewhat larger A). The Bogoliubov 
de Gennes excitation spectrum shown in Figs. [T47 a) and !14f b) is rich, with a series of crossings and avoided crossings. 
Here, we focus on the large D regime. Figures EJa) and[TJJb) show that the lowest k = 3 mode becomes soft first, 
followed by the lowest k = 2, k = 1 and k = modes. This indicates that the collapse is triggered by a mode with 
non- vanishing azimuthal quantum number, similarly to the case with vanishing barrier height [l5[. The D value, 
D w 163, at which the k = 3 mode becomes soft is about 16% smaller than the D value at which the stationary GP 
equation first supports negative energy solutions. This implies an appreciable reduction of the S>o islands shown in 
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FIG. 12: (Color online) Excitation spectrum obtained by solving the Bogoliubov de Gennes equations as a function of D for 
A = 12E Z , b — 0.2a z and A = 5. Solid and dashed lines in panel (a) show the real parts of the Bogoliubov de Gennes excitation 
frequencies for k = and 1, while solid and dashed lines in panel (b) show those for k = 2 and 3. 





FIG. 13: (Color online) Bogoliubov de Gennes eigenmodes f(p, z) corresponding to the lowest non-vanishing k = frequency 
for A = 5, A = 12E Z , b = 0.2a z and (a) D — 13.42 and (b) D — 21.87. The contours are chosen equidistant, with solid and 
dashed lines corresponding to positive and negative values of /. The dash-dotted lines indicate the nodal line of /. 
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FIG. 14: (Color online) Excitation spectrum obtained by solving the Bogoliubov de Gennes equations as a function of D for 
A = 12E Z , b — 0.2a z and A = 7. Solid and dashed lines in panel (a) show the real parts of the Bogoliubov de Gennes excitation 
frequencies for k = and 1, while solid and dashed lines in panel (b) show those for k = 2 and 3. 

Figs. OQ and 

It is worth emphasizing at this point that the fact that the k = 3 mode becomes soft first is not specific to the 
double-well geometry considered here; in fact, decay triggered by the k = 3 mode has also been found for pancake- 
shaped dipolar gases without barrier (see Ref. [l5[ and below) . While Fig. Q3] shows an example where the k = 3 
mode becomes soft first, we find that the collapse can — again, just as in the case of vanishing barrier [l5| — also be 
triggered by other finite k modes. For example, for a somewhat smaller barrier height but the same barrier width 
and aspect ratio as in Fig. [14] (i.e., for A — 9E Z , b = Q.2a z and A = 7), we find that the k = 2 mode becomes soft 
first, followed by the k = 3 and k = 1 modes (our calculations included modes k — through 4). In the following, we 
analyze the collapse triggered by the k = 3 mode in more detail. 

Figure [15] shows the eigenmodes associated with the two lowest k = 3 Bogoliubov de Gennes frequencies for 
D = 119.1 [Figs. HHa) and [T^b)] and D = 161.4 [Figs. [1%) andHgd)]. Figure [15] shows that the eigenmode 
corresponding to the lowest non-vanishing k — 3 frequency has no nodal line prior to the avoided crossing at D « 140 
[see Fig. 115( a)] but one nodal line, which can be parametrized roughly by p w 6a z (or, equivalently by p ss 2.6a p ), just 
prior to collapse [see Fig. [TBlc)]. Importantly, the eigenmode remains symmetric with respect to z = even close to 
collapse and exhibits two equivalent extrema at positive and negative z. This suggests that the collapse is associated 
with a total of six density peaks, three located on a ring of the red blood cell located in the left well and three located 
on a ring of the red blood cell located in the right well. The collapse can thus be characterized as local, with the local 
character arising from i) the angular roton like nature of the instability and ii) the presence of the Gaussian barrier. 

It is interesting to compare the eigenmodes for systems with vanishing and finite barrier in the regime where the 
first Bogoliubov de Gennes mode becomes soft. To this end, Fig. [TBI shows the eigenmodes corresponding to the two 
lowest non- vanishing k — 3 frequencies for two different D values, i.e., for D — 66.14 and D = 81.49, and A = 7 and 
A = 0. The latter D value corresponds to that investigated in Fig. 211 of Ref. [la ]. While we find that our excitation 
spectrum for A = (not shown here) agrees with Fig. 2IIb of Ref. [15| . our eigenmode corresponding to the lowest 
non- vanishing k — 3 frequency differs. In fact, we find that the eigenmode shown in Fig. 2IIc for D = 81.49 (rJlj 
corresponds to the second lowest and not to the lowest non-vanishing k — 3 Bogoliubov de Gennes eigen frequency 
as stated in Ref. [lj| (62J. The eigenmode shown in Fig. UM c) for A — and D — 81.49 has a nodal line very similar 
to that shown in Fig. [TSTc) for A = 12E Z and D = 161.4, i.e., for a D value that is roughly twice as large as that for 
A = 0. The main difference between the two eigenmodes is that the latter is characterized by extrema at z « ±a z as 
opposed to z = 0. 
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FIG. 15: (Color online) Bogoliubov de Gennes eigenmodes f(p,z) corresponding to the two lowest non-vanishing k — 3 
frequencies for A = 7, A — 12E Z and b = 0.2a z . Panels (a) and (b) show the eigenmodes corresponding to the lowest and 
second lowest non-vanishing k — 3 Bogoliubov de Gennes frequencies for D = 119.1 (i.e., prior to the avoided crossing at 
D « 140), while panels (c) and (d) show those for D = 161.4 (i.e., just prior to the lowest k — 3 mode becoming soft). The 
contours are chosen equidistant, with solid and dashed lines corresponding to positive and negative values of /. The dash-dotted 
lines indicate the nodal lines of /. 




FIG. 16: (Color online) Bogoliubov de Gennes eigenmodes f(p,z) corresponding to the two lowest non-vanishing k — 3 
frequencies for A = 7, A — and b — 0.2a z . Panels (a) and (b) show the eigenmodes corresponding to the lowest and second 
lowest non-vanishing k — 3 Bogoliubov de Gennes frequencies for D — 66.14, while panels (c) and (d) show those for D = 81.49 
(i.e., just prior to the lowest k = 3 mode becoming soft). The contours are chosen equidistant, with solid and dashed lines 
corresponding to positive and negative values of /. The dash-dotted lines indicate the nodal lines of /. 
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VI. SUMMARY AND OUTLOOK 

This paper presents a detailed mean-field analysis of a purely dipolar BEC in a cylindrically symmetric external 
confining potential with repulsive Gaussian barrier centered at z — 0. The dipoles are assumed to be aligned along 
one of the symmetry axes of the confining potential, which can be realized experimentally through the application of 
an external field. We have investigated the behaviors of the system as functions of the dimensionless parameter D, 
which is defined as the product of the number of particles and the square of the magnitude of the dipolc moment d, 
the aspect ratio A and, in a few selected cases, the height A of the Gaussian barrier. Throughout, the barrier width 
b and the s-wave scattering length a s were kept fixed at b — 0.2a z and a s = 0, respectively. The energetics and the 
density profiles obtained by solving the stationary GP equation have been discussed and the onset of the mechanical 
instability has been analyzed. Additional insights were gained from a dynamical stability analysis, which is based on 
time-evolving a given initial state or on the Bogoliubov de Gennes excitation spectrum. The latter was complemented 
by a detailed analysis of selected Bogoliubov de Gennes eigenmodes. 

For sufficiently small aspect ratios, i.e., for cigar-shaped harmonic traps in which the dipoles are aligned along 
the weak confinement direction, the energetically lowest lying stationary ground state solution is either symmetric 
or symmetry-broken. As in the case of s-wave interacting BECs, the appearance of these solutions can be explained 
qualitatively within a two-mode model. For a critical mean- field strength D cr , the system becomes unstable. The 
Bogoliubov de Gennes framework reveals that the collapse occurs globally, i.e., within a single well, with the decay 
being triggered by a mode with vanishing azimuthal quantum number. 

As the aspect ratio A increases, the region of the phase diagram where the energetically lowest lying stationary GP 
solution is symmetry-broken vanishes. For sufficiently large aspect ratio A, the symmetric solutions fall into one of 
two classes: The system's density is characterized by a density maximum at p = (this is referred to as So type) 
or by a density maximum at p > (this is referred to as S>o type). The latter class of solutions occupies a fairly 
small region of the phase diagram and only occurs near the dynamical instability line and only for certain aspect 
ratios A and barrier heights A. For a barrier height of A — 12E Z and A = 7, e.g., the solution is of type S>o prior 
to collapse. Correspondingly, we find that the collapse occurs locally through a mode with finite azimuthal quantum 
number (A; = 3), with density spikes emerging in six different regions of the trap. Three of these density peaks grow in 
the left well and three in the right well. Furthermore, each of the three density peaks lies on a ring that is associated 
with the density maximum of the ground state solution at p > 0. The instability can, as in the case of a vanishing 
Gaussian barrier, be characterized as an angular roton instability. We note, however, that the radial degrees of 
freedom also play a role, i.e., that the Bogoliubov de Gennes eigenmodes prior to collapse contain radial nodal lines. 
For pancake-shaped trapping geometries (i.e., for A > 1), this paper primarily explored the (D,X) parameter space 
for fixed A. We find that the double well system exhibits a number of novel and rich stability characteristics as the 
barrier height A is varied; these studies will be reported on in a forthcoming article [63l |. 

Our theoretical predictions for purely dipolar BECs presented in this paper can be tested experimentally by loading 
a dipolar BEC such as a Cr BEC into a double well potential and by tuning the s-wave scattering length to zero 
through the application of an external magnetic field in the vicinity of a magnetic Feshbach resonance. Following the 
spirit of the double- well experiments for s-wave interacting BECs (see, e.g., Ref. HU), ^ should be possible to study 
the transition from the Josephson tunneling regime to the macroscopic quantum self-trapping regime by loading the 
double well system with varying number of particles and varying population difference in the left and right wells. 
For pancake-shaped harmonic confinement, we suggest an experimental sequence that would allow the stability lines 
discussed in the context of Fig. [5] to be probed. We suggest to increase the radial trapping frequency uj p (and to thus 
decrease A) for fixed A and to monitor the loss of atoms from the trap. This scenario corresponds to approaching the 
instability line in Fig. [5] vertically from above. By repeating this experiment for condensates with varying number of 
particles, the different collapse mechanisms associated with the So regions and the S>o islands could be probed (see 
also Ref. (g|). 

In the future, it will be interesting to investigate how the behaviors of the system change as a function of the 
"spacing" between the left and the right well, i.e., as a function of the barrier widths b. The present study covers the 
regime where the "spacing" between the left and the right well is comparatively small, i.e., where it is comparable 
to the axial confinement length, and where the system is described by one macroscopic wave function. Our approach 
for pancake-shaped confinement with Gaussian barrier is distinctly different from other recent studies of multi-layer 
(quasi-)two-dimensional dipolar BECs HH, which assume that the dipoles in neighboring wells feel each other but 
that the distance between the neighboring wells is so large that each dipole can be assigned to a specific well. In this 
case, the system has been discretized, leading to a coupled set of equations that have been solved self-consistently. 
It will be interesting to extend the present study of pancake-shaped two-well systems to a regime where comparisons 
with a discretized description become meaningful. It will also be interesting to extend the present work to multi-well 
traps. In the regime of small aspect ratio, e.g., a three- or four-well system might lead to interesting dynamics that 
can be controlled by varying the onsite and the offsite interactions. In this case, a multi-mode analysis suggests itself 
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as a first starting point. 
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